home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / pascal / complxo.exe / COMPLEXO.PAS < prev    next >
Pascal/Delphi Source File  |  1992-01-15  |  23KB  |  729 lines
  1. {$N+,E+} {Use math coprocessor, if any, emulate otherwise.}
  2. UNIT ComplexOps;
  3.  
  4.  {This UNIT provides complex arithmetic and transcendental functions.
  5.  
  6.   (C) Copyright 1990, 1992, Earl F. Glynn, Overland Park, KS.  Compuserve 73257,3527.
  7.   All rights reserved.  This program may be freely distributed only for
  8.   non-commercial use.
  9.  
  10.   Some ideas in this UNIT were borrowed from "A Pascal Tool for Complex
  11.   Numbers", Journal of Pascal, Ada, & Modula-2, May/June 1985, pp. 23-29.
  12.   Many complex formulas were taken from Chapter 4, "Handbook of Mathematical
  13.   Functions" (Ninth Printing), Abramowitz and Stegun (editors), Dover, 1972.}
  14.  
  15. INTERFACE
  16.  
  17.   TYPE
  18.     RealType    = DOUBLE;
  19.     ComplexForm = (polar,rectangular);
  20.     Complex =
  21.       RECORD
  22.         CASE form:  ComplexForm OF
  23.           rectangular:  (x,y    :  RealType);  {z = x + y*i}
  24.           polar      :  (r,theta:  RealType);  {z = r*CIS(theta)}
  25.       END;                 {where CIS(theta) = COS(theta) + SIN(theta)*i}
  26.                            {      theta = -PI..PI (in canonical form)}
  27.  
  28.   CONST
  29.     MaxTerm    : BYTE     = 35;
  30.     EpsilonSqr : RealType = 1.0E-20;
  31.     Infinity   : RealType = 1.0E25;    {virtual infinity}
  32.  
  33.                              {complex number definition/conversion/output}
  34.   PROCEDURE CConvert (VAR z:  Complex; f:  ComplexForm);
  35.   PROCEDURE CSet (VAR z:  Complex; a,b:  RealType; f:  ComplexForm);
  36.   FUNCTION  CStr (z:  Complex; w,d:  BYTE; f:  ComplexForm):  STRING;
  37.  
  38.                              {complex arithmetic}
  39.   PROCEDURE CAdd  (VAR z:  Complex; a,b:  Complex);    {z = a + b}
  40.   PROCEDURE CDiv  (VAR z:  Complex; a,b:  Complex);    {z = a / b}
  41.   PROCEDURE CMult (VAR z:  Complex; a,b:  Complex);    {z = a * b}
  42.   PROCEDURE CSub  (VAR z:  Complex; a,b:  Complex);    {z = a - b}
  43.   PROCEDURE CNeg  (VAR z:  Complex; a  :  Complex);    {z = -a   }
  44.  
  45.                              {complex natural log, exponential}
  46.   PROCEDURE CLn   (VAR fn :  Complex; z:  Complex);   {fn  = ln(z) }
  47.   PROCEDURE CExp  (VAR z  :  Complex; a:  Complex);   {z   = exp(a)}
  48.   PROCEDURE CPwr  (VAR z  :  Complex; a,b:  Complex); {z   = a^b   }
  49.  
  50.                              {complex trig functions}
  51.   PROCEDURE CCos  (VAR z:  Complex; a:  Complex);     {z   = cos(a)}
  52.   PROCEDURE CSin  (VAR z:  Complex; a:  Complex);     {z   = sin(a)}
  53.   PROCEDURE CTan  (VAR z:  Complex; a:  Complex);     {z   = tan(a)}
  54.  
  55.   PROCEDURE CSec  (VAR z:  Complex; a:  Complex);     {z   = sec(a)}
  56.   PROCEDURE CCsc  (VAR z:  Complex; a:  Complex);     {z   = csc(a)}
  57.   PROCEDURE CCot  (VAR z:  Complex; a:  Complex);     {z   = cot(a)}
  58.  
  59.                              {complex hyperbolic functions}
  60.   PROCEDURE CCosh (VAR z:  Complex; a:  Complex);     {z   = cosh(a)}
  61.   PROCEDURE CSinh (VAR z:  Complex; a:  Complex);     {z   = sinh(a)}
  62.   PROCEDURE CTanh (VAR z:  Complex; a:  Complex);     {z   = tanh(a)}
  63.  
  64.   PROCEDURE CSech (VAR z:  Complex; a:  Complex);     {z   = sech(a)}
  65.   PROCEDURE CCsch (VAR z:  Complex; a:  Complex);     {z   = csch(a)}
  66.   PROCEDURE CCoth (VAR z:  Complex; a:  Complex);     {z   = coth(a)}
  67.  
  68.                              {miscellaneous complex functions}
  69.   FUNCTION  CAbs (z:  Complex):  RealType;                 {CAbs = |z|}
  70.   FUNCTION  CAbsSqr (z:  Complex):  RealType;           {CAbsSqr = |z|^2}
  71.   PROCEDURE CIntPwr (VAR z:  Complex; a:  Complex; n:  INTEGER); {z = a^n}
  72.   PROCEDURE CRealPwr (VAR z:  Complex; a:  Complex; x:  RealType); {z = a^x}
  73.   PROCEDURE CConjugate (VAR z:  Complex; a:  Complex);        {z = a*}
  74.   PROCEDURE CSqrt (VAR z:  Complex; a:  Complex);      {z = SQRT(a)}
  75.   PROCEDURE CRoot (VAR z:  Complex; a:  Complex; k,n:  WORD); {z = a^(1/n)}
  76.  
  77.                              {complex Bessel functions of order zero}
  78.   PROCEDURE CI0   (VAR sum:  Complex; z:  Complex);  {sum = I0(z)}
  79.   PROCEDURE CJ0   (VAR sum:  Complex; z:  Complex);  {sum = J0(z)}
  80.  
  81.   PROCEDURE CLnGamma (VAR z:  Complex; a:  Complex);
  82.   PROCEDURE CGamma   (VAR z:  Complex; a:  Complex);
  83.  
  84.                                   {treat "fuzz" of real numbers}
  85.   PROCEDURE CDeFuzz (VAR z:  Complex);
  86.   FUNCTION  DeFuzz (x:  RealType):  RealType;
  87.   PROCEDURE SetFuzz (value:  RealType);
  88.  
  89.                                   {miscellaneous}
  90.   FUNCTION FixAngle (theta:  RealType):  RealType;    {-PI < theta <= PI}
  91.   FUNCTION Pwr (x,y:  RealType):  RealType;    {Pwr = x^y}
  92.   FUNCTION Log10 (x:  RealType):  RealType;
  93.   FUNCTION LongMod (l1,l2:  LongInt):  LongInt;
  94.   FUNCTION Cosh (x:  RealType):  RealType;
  95.   FUNCTION Sinh (x:  RealType):  RealType;
  96.  
  97. IMPLEMENTATION
  98.  
  99.   VAR
  100.     fuzz     :  RealType;
  101.     Cone     :  Complex;
  102.     Cinfinity:  Complex;
  103.     Czero    :  Complex;
  104.     hLnPI    :  RealType;
  105.     hLn2PI   :  RealType;
  106.     ln2      :  RealType;
  107.  
  108.                              {complex number definition/conversion/output}
  109.   PROCEDURE CConvert (VAR z:  Complex; f:  ComplexForm);
  110.     VAR a:  Complex;
  111.   BEGIN
  112.     IF   z.form = f
  113.     THEN CDeFuzz (z)
  114.     ELSE BEGIN
  115.       CASE z.form OF
  116.         polar:            {polar-to-rectangular conversion}
  117.           BEGIN
  118.             a.form := rectangular;
  119.             a.x := z.r * COS(z.theta);
  120.             a.y := z.r * SIN(z.theta)
  121.           END;
  122.         rectangular:      {rectangular-to-polar conversion}
  123.           BEGIN
  124.             a.form := polar;
  125.             IF   DeFuzz(z.x) = 0.0
  126.             THEN BEGIN
  127.               IF   DeFuzz(z.y) = 0.0
  128.               THEN BEGIN
  129.                 a.r     := 0.0;
  130.                 a.theta := 0.0
  131.               END
  132.               ELSE
  133.                 IF   z.y > 0.0
  134.                 THEN BEGIN
  135.                   a.r := z.y;
  136.                   a.theta := 0.5*PI
  137.                 END
  138.                 ELSE BEGIN
  139.                   a.r := -z.y;
  140.                   a.theta := -0.5*PI
  141.                 END
  142.             END
  143.             ELSE BEGIN
  144.               a.r := CAbs(z);
  145.               a.theta := ARCTAN(z.y/z.x);   {4th/1st quadrant -PI/2..PI/2}
  146.               IF   z.x < 0.0                {2nd/3rd quadrants}
  147.               THEN
  148.                 IF   z.y >= 0.0
  149.                 THEN a.theta :=  PI + a.theta {2nd quadrant:  PI/2..PI}
  150.                 ELSE a.theta := -PI + a.theta {3rd quadrant: -PI..-PI/2}
  151.             END
  152.           END;
  153.       END;
  154.       CDeFuzz (a);
  155.       z := a
  156.     END
  157.   END {CConvert};
  158.  
  159.   PROCEDURE CSet (VAR z:  Complex; a,b:  RealType; f:  ComplexForm);
  160.   BEGIN
  161.     z.form := f;
  162.     CASE f OF
  163.       polar:
  164.         BEGIN
  165.           z.r := a;
  166.           z.theta := b
  167.         END;
  168.       rectangular:
  169.         BEGIN
  170.           z.x := a;
  171.           z.y := b
  172.         END;
  173.     END
  174.   END {CSet};
  175.  
  176.   FUNCTION  CStr (z:  Complex; w,d:  BYTE; f:  ComplexForm):  STRING;
  177.     VAR s1,s2:  STRING;
  178.   BEGIN
  179.     CConvert (z,f);
  180.     CASE f OF
  181.       polar:
  182.         BEGIN
  183.           Str (z.r:w:d, s1);
  184.           Str (z.theta:w:d, s2);
  185.           CStr := s1+'*CIS('+s2+')'
  186.         END;
  187.       rectangular:
  188.         BEGIN
  189.           Str (z.x:w:d, s1);
  190.           Str (ABS(z.y):w:d, s2);
  191.           IF   z.y >= 0
  192.           THEN CStr := s1+' +'+s2+'i'
  193.           ELSE CStr := s1+' -'+s2+'i'
  194.         END
  195.     END
  196.   END {CStr};
  197.  
  198.                                   {complex arithmetic}
  199.   PROCEDURE CAdd  (VAR z:  Complex; a,b:  Complex);    {z = a + b}
  200.   BEGIN                                               {complex addition}
  201.     CConvert (a,rectangular);
  202.     CConvert (b,rectangular);
  203.     z.form := rectangular;
  204.     z.x := a.x + b.x;   {real part}
  205.     z.y := a.y + b.y;   {imaginary part}
  206.   END {CAdd};
  207.  
  208.   PROCEDURE CDiv  (VAR z:  Complex; a,b:  Complex);    {z = a / b}
  209.     VAR temp:  RealType;
  210.   BEGIN
  211.     CConvert (b,a.form);    {arbitrarily convert one to type of other}
  212.     z.form := a.form;
  213.     CASE a.form OF
  214.       polar:
  215.         BEGIN
  216.           z.r := a.r / b.r;
  217.           z.theta := FixAngle(a.theta - b.theta)
  218.         END;
  219.       rectangular:
  220.         BEGIN
  221.           temp := SQR(b.x) + SQR(b.y);
  222.           z.x := (a.x*b.x + a.y*b.y) / temp;
  223.           z.y := (a.y*b.x - a.x*b.y) / temp
  224.         END
  225.     END
  226.   END {CDiv};
  227.  
  228.   PROCEDURE CMult (